/*
 * CannonicalSystem.cpp
 *
 *  Created on: Jan 9, 2012
 *      Author: allan
 */

#include "CannonicalSystem.h"

// Eigen includes
#include <Eigen/LU>

CannonicalSystem::CannonicalSystem()
{}

CannonicalSystem::CannonicalSystem(const Multiphase& multiphase,
	const ReactionSystem& equilibrium_reactions,
	const vector<string>& primary_species,
	const vector<string>& equilibrium_species)
: multiphase(multiphase), equilibrium_reactions(equilibrium_reactions),
  primary_species(primary_species), equilibrium_species(equilibrium_species)
{
	indexes_primary_species = multiphase[primary_species];
	indexes_equilibrium_species = multiphase[equilibrium_species];
	
	uej = equilibrium_reactions.AssembleStoichiometricMatrix(primary_species);
	uee = equilibrium_reactions.AssembleStoichiometricMatrix(equilibrium_species);
	
	MatrixXd uee_inv = uee.lu().inverse();
	
	vej = - uee_inv * uej;
}

CannonicalSystem::CannonicalSystem(const Multiphase& multiphase,
	const ReactionSystem& equilibrium_reactions,
	const ReactionSystem& kinetic_reactions,
	const vector<string>& primary_species,
	const vector<string>& equilibrium_species,
	const vector<string>& kinetic_species)
: multiphase(multiphase), equilibrium_reactions(equilibrium_reactions),
  kinetic_reactions(kinetic_reactions), primary_species(primary_species),
  equilibrium_species(equilibrium_species), kinetic_species(kinetic_species)
{
	indexes_primary_species = multiphase[primary_species];
	indexes_equilibrium_species = multiphase[equilibrium_species];
	indexes_kinetic_species = multiphase[kinetic_species];
	
	uej = equilibrium_reactions.AssembleStoichiometricMatrix(primary_species);
	uee = equilibrium_reactions.AssembleStoichiometricMatrix(equilibrium_species);
	
	ukj = kinetic_reactions.AssembleStoichiometricMatrix(primary_species);
	ukk = kinetic_reactions.AssembleStoichiometricMatrix(kinetic_species);
	
	MatrixXd uee_inv = uee.lu().inverse();
	MatrixXd ukk_inv = ukk.lu().inverse();
	
	vej = - uee_inv * uej;
	vkj = - ukk_inv * ukj;
}

const unsigned
CannonicalSystem::GetNumSpecies() const
{
	return primary_species.size() + equilibrium_species.size() + kinetic_species.size();
}

const unsigned
CannonicalSystem::GetNumPrimarySpecies() const
{
	return primary_species.size();
}

const unsigned
CannonicalSystem::GetNumEquilibriumSpecies() const
{
	return equilibrium_species.size();
}

const unsigned
CannonicalSystem::GetNumKineticSpecies() const
{
	return kinetic_species.size();
}

const Multiphase&
CannonicalSystem::GetMultiphase() const
{
	return multiphase;
}

const ReactionSystem&
CannonicalSystem::GetEquilibriumReactions() const
{
	return equilibrium_reactions;
}

const ReactionSystem&
CannonicalSystem::GetKineticReactions() const
{
	return kinetic_reactions;
}

const vector<string>&
CannonicalSystem::GetPrimarySpecies() const
{
	return primary_species;
}

const vector<string>&
CannonicalSystem::GetEquilibriumSpecies() const
{
	return equilibrium_species;
}

const vector<string>&
CannonicalSystem::GetKineticSpecies() const
{
	return kinetic_species;
}

const vector<Index>&
CannonicalSystem::GetIndexesPrimarySpecies() const
{
	return indexes_primary_species;
}

const vector<Index>&
CannonicalSystem::GetIndexesEquilibriumSpecies() const
{
	return indexes_equilibrium_species;
}

const vector<Index>&
CannonicalSystem::GetIndexesKineticSpecies() const
{
	return indexes_kinetic_species;
}

void
CannonicalSystem::SetCompositionPrimarySpecies(const VectorXd& nj, VectorXd& n) const
{
	const unsigned Nj = GetNumPrimarySpecies();
	
	for(unsigned j = 0; j < Nj; ++j)
		n[indexes_primary_species[j]] = nj[j];
}

void
CannonicalSystem::SetCompositionEquilibriumSpecies(const VectorXd& ne, VectorXd& n) const
{
	const unsigned Ne = GetNumEquilibriumSpecies();
	
	for(unsigned e = 0; e < Ne; ++e)
		n[indexes_equilibrium_species[e]] = ne[e];
}

void
CannonicalSystem::SetCompositionKineticSpecies(const VectorXd& nk, VectorXd& n) const
{
	const unsigned Nk = GetNumKineticSpecies();
	
	for(unsigned k = 0; k < Nk; ++k)
		n[indexes_kinetic_species[k]] = nk[k];
}

const MatrixXd
CannonicalSystem::StoichiometricMatrixUej() const
{
	return uej;
}

const MatrixXd
CannonicalSystem::StoichiometricMatrixUee() const
{
	return uee;
}

const MatrixXd
CannonicalSystem::StoichiometricMatrixUkj() const
{
	return ukj;
}

const MatrixXd
CannonicalSystem::StoichiometricMatrixUkk() const
{
	return ukk;
}

const MatrixXd
CannonicalSystem::StoichiometricMatrixVej() const
{
	return vej;
}

const VectorXd
CannonicalSystem::CompositionPrimarySpecies(const VectorXd& n) const
{
	const unsigned Nj = GetNumPrimarySpecies();
	
	VectorXd nj(Nj);
	
	for(unsigned j = 0; j < Nj; ++j)
		nj[j] = n[indexes_primary_species[j]];
	
	return nj;
}

const VectorXd
CannonicalSystem::CompositionEquilibriumSpecies(const VectorXd& n) const
{
	const unsigned Ne = GetNumEquilibriumSpecies();
	
	VectorXd ne(Ne);
	
	for(unsigned e = 0; e < Ne; ++e)
		ne[e] = n[indexes_equilibrium_species[e]];
	
	return ne;
}

const VectorXd
CannonicalSystem::CompositionKineticSpecies(const VectorXd& n) const
{
	const unsigned Nk = GetNumKineticSpecies();
	
	VectorXd nk(Nk);
	
	for(unsigned k = 0; k < Nk; ++k)
		nk[k] = n[indexes_kinetic_species[k]];
	
	return nk;
}

const VectorXd
CannonicalSystem::KineticRatesPrimarySpecies(const VectorXd& rk) const
{
	return ukj.transpose() * rk;
}

const VectorXd
CannonicalSystem::KineticRatesKineticSpecies(const VectorXd& rk) const
{
	return ukk.transpose() * rk;
}

const VectorXd
CannonicalSystem::TotalCompositionPrimarySpecies(const VectorXd& n) const
{
	const VectorXd nj = CompositionPrimarySpecies(n);
	const VectorXd ne = CompositionEquilibriumSpecies(n);
	
	return nj + vej.transpose() * ne;
}

